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Abstract 

Irregular structure in planetary rings is often attributed to the intrinsic instabilities 
of a homogeneous state undergoing Keplerian shear. Previously these have been 
analysed with simple hydrodynamic models. We instead employ a kinetic theory, 
in which we solve the linearised moment equations derived in Shu and Stewart 
1985 for a dilute ring. This facilitates an examination of velocity anisotropy and 
non-Newtonian stress, and their effects on the viscous and viscous/gravitational 
instabilities thought to occur in Saturn's rings. Because we adopt a dilute gas model, 
the applicability of our results to the actual dense rings of Saturn are significantly 
curtailled. Nevertheless this study is a necessary preliminary before an attack on 
the difficult problem of dense ring dynamics. 

We find the Shu and Stewart formalism admits analytic stability criteria for the 
viscous overstability, viscous instability, and thermal instability. These criteria are 
compared with those of a hydrodynamic model incorporating the effective viscos- 
ity and cooling function computed from the kinetic steady state. We find the two 
agree in the 'hydrodynamic limit' (i.e. many collisions per orbit) but disagree when 
collisions are less frequent, when we expect the viscous stress to be increasingly 
non-Newtonian and the velocity distribution increasingly anisotropic. In particular, 
hydrodynamics predicts viscous overstability for a larger portion of parameter space. 
We also numerically solve the linearised equations of the more accurate Goldreich 
and Tremaine 1978 kinetic model and discover its linear stability to be qualitatively 
the same as that of Shu and Stewart's. Thus the simple collision operator adopted 
in the latter would appear to be an adequate approximation for dilute rings, at least 
in the linear regime. 
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1 Introduction 



This paper examines the coUisional dynamics of differentially rotating, partic- 
ulate disks, and addresses specifically the intrinsic mechanisms which might 
induce the generation of irregular, fine-scale structure, such as that exhibited 
by Saturn's B-ring. Voyager, and Cassini more recently, revealed that these 
essentially axisymmetric variations appear on a very wide range of length- 
scale- from the hmit of resolution (~ 100m) to several hundred kilometres 
(Horn and Cuzzi 1995, Porco et al. 2005 ). Though some proportion of B-ring 
structure may be caused by non-dynamical effects such as variations in albedo 
and phase function- see Cuzzi and Estrada 1996, for example) and from grav- 
itational interactions with neighbouring satellites (Thiessenhusen et al. 1995), 
it is believed that the majority corresponds to changes in optical thickness 
resulting from a collective dynamics (Tremaine 2003). 

This work studies dilute particulate systems, and as such may not be directly 
applicable to the saturnian rings. It is nevertheless an informative first stop 
before we embark on the challenging stability analysis of a dense granular 
gas. Other than presenting a general framework in which such a calculation 
can be tackled, the dilute gas formalism permits us to isolate analytically the 
interesting effects of anisotropy and non-Newtonian viscous stress upon the 
local axisymmetric instabilities thought to be manifest in planetary rings. 

Since Voyager first reported the radial stratification of Saturn's rings, theo- 
reticians have advocated a number of possible causes. These have inchidcd: 
ballistic transport (Duriscn 1995), electromagnetic effects (Goertz and Morfill 
1988), and interleaved shearing and shear-free zones (Tremaine 2003). Others 
have looked to the local instabiUties of viscous fluid disks, interpreting the 
disordered state observed as the saturated endpoint of their nonlinear evolu- 
tion. The 'viscous instability' was the first proposed (Lin and Bodenheimer 
1981, Ward 1981, and Lukkari 1981). Essentially a monotonic 'clumping', it 
is associated with an outward angular momentum fiux which decreases with 
surface density: diva) /da < (where v is kinematic viscosity and a surface 
mass density). Though a dilute ring's viscosity depends on optical depth in 
a manner which promises the existence of such an instability (Goldreich and 
Tremaine 1978, Shu and Stewart 1985), Saturn's rings are most likely 'dense', 
and theoretical and numerical N-body studies have revealed that such rings 
do not manifest the appropriate viscosity for this instability to develop (Araki 
and Tremaine 1986, and Wisdom and Tremaine 1988). 

The 'viscous overstability' was first examined in the context of accretion disks 
by Kato (1978), and, as the name suggests, originates in overcompensation 
by the system's restoring forces: the stress oscillation which accompanies the 
epicyclic response in an acoustic-inertial wave will force the system back to 
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equilibrium so strongly that it will 'overshoot'. The mechanism relics on: a) 
the synchronisation of the viscous stress's oscillations with those of density, 
and b) the viscous stress increasing sufficiently in the compressed phase. In 
hydrodynamics only the latter consideration is relevant, which furnishes the 
criterion for overstability: (5 = {din 1/ / din a) > (3*, where /3* is a number 
dependent on the thermal properties of the ring (Schmit and Tscharnuter 
1995). 

The viscous overstability has been a favoured explanation for smaller scale 
B-ring structure in recent years, a status stemming, primarily, from the vis- 
cosity profiles computed by Araki and Tremaine's 1986 dense gas model and 
Wisdom and Tremaine's 1988 particle simulations. Both appear to satisfy the 
above criterion. Consequently the linear behaviour of the instability has been 
thoroughly examined, though only within a hydrodynamic framework (Schmit 
and Tscharnuter 1995, Spahn et. al. 2000, and Schmidt et al. 2001). In ad- 
dition, Schmidt and Salo (2003) have constructed a weakly nonfinear theory, 
and the overstability' s long-term, nonlinear behaviour has been numerically 
studied by Schmit and Tscharnuter (1999). An isothermal model was adopted 
in both cases. 

However, fluid mechanics is only one of several theoretical approaches thathave 
been deployed to capture the behaviour of planetary rings: N-body simula- 
tions, generalisations of stellar dynamics, and kinetic theory. This plurality is 
surely indicative of how the field of planetary ring dynamics falls uncomfort- 
ably between the more familiar frameworks of classical physics. Our analysis 
will predominantly focus upon the kinetic approach, and solve the set of mo- 
ment equations derived from the kinetic theories proposed in Shu and Stewart 
1985 (hereafter referred to as 'SS85' ) and in Goldreich and Tremaine 1978 
(referred to as 'GT78'). These formafisms have hitherto been used mainly to 
establish equilibrium solutions. This paper goes the next step by exploring 
their linear theory. 

Such an analysis is more involved than the analogous hydrodynamic calcula- 
tion; however, the latter's adoption of the Navier-Stokes stress model intro- 
duces two assumptions which may be inappropriate in the ring context and 
whose consequences are instructive to investigate. Firstly, the Navier-Stokes 
model presumes the distribution of the particles' velocity dispersion to be 
nearly isotropic. In the regime of many collisions per orbit this is an acceptable 
supposition, as collisions scatter particles randomly on the average. However 
if the collision rate, cUc, is of the same order as the orbital frequency, fl, (as it 
is presumed to be in Saturn's rings) this need not be true. Secondly, hydrody- 
namics assumes an 'instantaneous' (local in time) relationship between stress 
and strain. This may not hold when ujc ~ fl. Generally the viscous stress pos- 
sesses a relaxation time of order l/cUc, which in this regime will be comparable 
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to the dynamic time scale. Thus the immediate history of the stress cannot be 
ignored. Its inchision should have the most impact on the stability of oscillat- 
ing modes, especially the overst ability, it depending on the synchronisation of 
the stress and density oscillations. A kinetic model can address both issues, 
accounting for anisotropy within an appropriate collision term and providing 
a straightforward way, by the taking of moments, to generate dynamical equa- 
tions for the viscous stress. Another advantage is that a kinetic model lets us 
explicitly include the microphysics of particle-particle interactions and thence 
to potentially model a larger set of the physical mechanisms at play (such as 
collisions, irregular surfaces, spin , size distribution, etc). It also narrows the 
scope of our simplifying assumptions to the particulars of collisions between 
spheres of ice, which have been observed in the laboratory (see Bridges et al. 
1984, Hatzes et al. 1988, amongst others) . 

The paper is organised as follows. In Section 2 we present a brief derivation 
of the set of closed, vertically averaged moment equations in a local model, 
the shearing sheet, for a general kinetic equation. We then specialise to the 
SS85 kinetic model and briefly exhibit its equilibrium characteristics. Section 
3 contains the linear theory of a general kinetic equation and proves results 
pertaining to self-gravity and the thermal modes when in the limit of length- 
scales much longer than the disk height. We then derive stability criteria for the 
viscous instability and overstability specific to the SS85 model. Comparisons 
with an analogous calculation for a hydrodynamic model and the GT78 model 
are performed in Sections 4 and 5. In the former we isolate the effects of non- 
Newtonian stress on the linear theory, and in the latter observe how the SS85 
collision term approximates satisfactorily that of the more accurate triaxial 
Gaussian model. We draw our conclusions in Section 6. 



2 Kinetic Model 

2.1 Simplifying Assumptions 

The formulation of a suitable kinetic theory poses a number of difficulties, 
namely the closure of the moment hierarchy and the simpliffcation of the col- 
lision term. But the various approximations these require by no means cripple 
the approach. The most fundamental assumption we make is that our ring 
is composed exclusively of hard, identical, and indestructible spheres. In ad- 
dition we presume that the particles are non-spinning. These simplifications 
render the mathematics convenient (but see Stewart et al. 1984 and Shu and 
Stewart 1985 for justifications)- now the inelasticity of the collisions is quan- 
tified solely by the normal coefficient of restitution, £, which is the ratio of the 
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normal relative speed after and before a collision. It is generally a function of 
normal impact velocity, w„, and possibly other parameters like ambient tem- 
perature and pressure, particle size and mass (Hatzes et al. 1988, Dilley 1993) 

A number of laboratory experiments have been conducted in order to ascertain 
the collisional properties of ice spheres in saturnian conditions, in particular 
the relationship between the coefficient of restitution and Vn (see Bridges et al. 
1984, Hatzes et al. 1988, Supulver et al. 1995, and Dilley and Crawford 1996). 
Excepting Dilley, these authors successfully fit a step-wise power law to their 
data for coUisions sufficiently gentle and/or surfaces sufficiently frosted: 



where Vc and p are parameters contingent on the material properties of the 
ice balls and their ambient environment. Bridges' data admit p = 0.234 and 
Vc = 0.0077 cm s^^ for frosted particles of radius 2.5cm at a temperature of 
21 OK (significantly higher than the appropriate conditions). Hatzes finds p ~ 
0.20 and Vc — 0.025 cm s~^ for the case of frosted particles at 123K (however, 
for the case of smoother particles an exponential law provides a better fit). 
At slightly lower temperatures (~100K) Supulver et al. obtain p = 0.19 and 
Vc = 0.029 cm s~^ with their torsional pendulum fixed. Though the functional 
relationship (1) is useful, the coefficient of restitution varies considerably as the 
physical condition of the contact surface is more or less frosty or sublimated, as 
Hatzes et al. explore. Also the neglect of spin, the effects of glancing collisions, 
irregularly shaped surfaces, coagulation and erosion, mass transfer and the 
role of amorphous ice, further emphasises the simplicity of the collision model 
adopted. Studies of those processes we omit are found in Salo 1987, Araki 1991 
and Morishima and Salo 2006 (theory and simulations), and in McDonald et 
al. 1989, Hatzes et al. 1991, Supulver et al. 1995 and Supulver et al. 1997 
(experiments) . 

From this point on we denote by e the averaged coefficient of restitution, 
which is a function of the macroscopic variable, c, the velocity dispersion. For 
a general step-wise power law and a Maxwellian velocity distribution, it is 
straightforward to obtain an analytic expression for e averaged over collisions. 
Its functional dependence on c is, unsurprisingly, the 'smoothed' analogue of 
Eq.(l)'s piecewise dependence on Vn. 

The rings we study are dilute. We realise this curtails the applicability of 
our results considerably as Saturn's rings, even the optically thin regions, 
are thought to exhibit important dense effects. Araki and Tremaine (1986) 
revealed significant collisional contributions to the kinetic equilibrium when 
collisions are sufficiently inelastic; and the experimentally derived elasticity 
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laws of (1) predict collisions are dissipative enough for this 'cold' regime to be 
widespread in Saturn's rings, given appropriate particle sizes (see Marouf et al. 
1983 and Zebker et al 1985). This conclusion has been confirmed by particle 
simulations, Wisdom and Tremaine 1988 and Salo 1991 the most notable. 
Moreover, both simulations and theoretical models show that the dynamics of 
small particles (for which nonlocal effects are less important) strongly couple to 
the dynamics of the largest particles (for which nonlocal effects are important), 
both in two size systems and in polydisperse disks exhibiting power-law size 
distributions akin to Saturn's (see Stewart et al. 1984, Salo 1987, Lukkari 
1989, Salo 1991, and Salo 1992b). 

Physically then, a dilute model best describes planetary rings in which the 
radii of the largest particles are of the order of centimetres, or tens of centime- 
tres. Estimates on the maximum size of particle radius are easy to establish 
from the dilute theory, though these are only necessary bounds: sufficient 
conditions for the applicability of the dilute model can only be determined 
from the dense theory. At lower optical depths, however, we expect the dilute 
estimates to serve as a reasonable guide. We obtain these by combining the 
famous equilibrium relation between optical thickness, r, and £, computed in 
GT78 (cf. Fig 7), with the elasticity laws embodied in (1). This determines 
the dependence of velocity dispersion on r. We then suppose that non-local 
effects may be neglected if a particle's r.m.s. speed, a/3c, is an order larger 
than the shear velocity across two particle radii, 3Qa, where a is particle ra- 
dius (Araki 1991)- the ratio of these quantities is a convenient measure of 
'nonlocality', but only for optical thicknesses less than about 1 or 2 (see Salo 
1991). Thereupon we compute an upper limit on a below which the dilute 
model is acceptable. If r = 0.2 at C-ring distances, flmax ranges from 31cm 
(for Supulver et al.^s data) to 5.3cm (for Bridges et a/.'s data). At A-ring 
distances for r = 0.5 the maximum particle radius is 33cm (Supulver et al.) 
or 6.3cm (Bridges et al). At B-ring distances and r — 1, the Bridges et al. 
data imply amax ~ 2.2cm, or, alternatively, Vc ~ 0.35 cm s^^ for metre-sized 
particles. This last estimate appears consistent with those simulations of Salo 
(1991) in which Vc was increased to 20 and 40 times vb = 0.01 cm/s (cf. his 
Fig.'s 5a, 5c and 5d). It is interesting to note the sensitivity of the critical 
a to the elasticity law adopted: relations which predict slightly less dissipa- 
tive impacts allow a substantially higher velocity dispersion, and subsequently 

^max- 

Our dilute ring is self-gravitating. Disk self-gravity may be decomposed into 
three effects. Firstly, inter-particle gravitational forces will lead to an addi- 
tional source of (elastic) scattering via gravitational encounters. They will 
also enhance the collision frequency of physical collisions (gravitational focus- 
ing). For either process to be appreciable, the velocity dispersion must be of 
order the particles' escape velocity, i.e. Gpa^ ~ c^, where G is the gravitation 
constant and p the mass density of the particle material. This can be recast 
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as {afl)/c ~ 1, for typical values of the saturnian system. As non-local effects 
are neglected, the left side of this scaling will be ^ 1, and hence gravitational 
scattering and focusing may be omitted. Secondly, the aggregate of all the 
particle attractions will increase the disk's vertical restoring force and further 
flatten the disk. Subsequently the collision frequency will increase, an effect 
which we include. Thirdly, there will arise a dynamical contribution to the 
disk's stability which issues from the extra term in the momentum equation. 
Though this term will be unimportant at equilibrium, being dominated by 
the planet's gravitational potential, it will be significant in the linear theory. 
As is well known, an inviscid fiuid disk collapses under the action of self- 
gravity for a Toomre parameter Q less than 1, its unstable modes growing 
on a band of intermediate wavelengths, with rotation and pressure stabilis- 
ing the long and short scales respectively (Toomre 1964, Julian and Toomre 
1966). A viscous disk however paints a slightly more complex picture. Then 
it is better to think of self-gravity as 'extending' the viscous instabilities we 
discuss into larger areas of parameter space. In these 'extensions' the instabil- 
ities grow only in a certain confined range of intermediate wavelength, unlike 
the non-self-gravitating case in which the longest wavelengths are the first to 
become viscously unstable and overstable (Schmit and Tscharnuter 1995). We 
shall examine this effect in some detail particularly for the viscous instability. 
We should mention also that self gravity induces transient, non-axisymmetric 
wakes, the analysis of which we omit (but see Salo 1992a). 

2.2 The Governing Moment Equations 

Consider a dilute gas of identical, smooth, inelastic spheres of mass m with 
phase space distribution /(x, v, t). The number of particles located in the 
volume (ix centred at x with velocities in the range dv centred at v at time 
t is defined as /(x, v, t)(ix(iv. Taking moments of / allows us to calculate 
the familiar macroscopic characteristics of the gas. Number density, n, bulk 
velocity, u, and velocity dispersion, c, are defined through 



The phase space distribution satisfies a kinetic equation which is distinguished 
by a collision operator, [df /dt)c, whose precise form refiects the coUisional 
microphysics. 

We situate our gas in a shearing sheet. This is a convenient representation of a 
differentially rotating disk in which a small patch, centred on a point moving 
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Fig. 1. The shearing sheet. The differential rotation is locally represented as a ro- 
tation superimposed upon a linear shear flow. 

on a circular orbit at r = tq, is represented as a sheet in uniform rotation, 
^{f^o)^z, and subjected to a linear shear flow, uq = —2Aoxey. The local 
rectilinear coordinates x and y point in the radial and azimuthal directions 
respectively, and = —^ro{d^l/dr)o; see Fig. 1. For the dynamical analysis 
we consider the sheet horizontally unbounded, at least compared to the vertical 
length scale of the disk. The model hence provides an excellent approximation 
for rings whose thickness is very small, such as Saturn's. 

In terms of the peculiar velocity, w = v — u, the kinetic equation for such a 
patch of ring is 



df , ,df 



dui dui 



9f _ (df\ ... 



where the force/unit mass is 



^ dx- ^eijkilj{wk + Uk). 

The appropriate centrifugal-gravitational potential of the planet is denoted 
by 0p, and the disk's gravitational potential by 0/). The tensor eijk is the 
alternating tensor and the angular velocity is O = r2oe^, where VLq = fl{ro). 
By multiplying (5) by 1, Wj and WiWj and then integrating over all w we derive 
the continuity equation, 

dtn + dk{nuk) = 0, (6) 

the equation of motion, 

n {dtUi + UkdkUi) = -2neijk^jUk - ndi{(f)p + (f)D) - djPij, (7) 

and the pressure tensor equation, 

dtPij + UkdkPij = -PikdkUj - PjkdkUi - PijdkUk - 2eikiflkPij 

— 2ejki^kPH — dkPijk + Qij, (8) 
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where pij and pijk are the second and third-order moments, 

Pij... = j WiWj . . . fdw, 
and the coUisional change in second moment is 

which is a yet to be specified function of the other field variables. For notational 
brevity we have set di — d/dxi. 




As we are interested in scales ^ H there is little point in solving for a disk's 
detailed vertical structure. We hence vertically integrate the moment equa- 
tions, a process which can be performed without much grief if it is assumed 
that Ux and Uy are effectively independent of z. The third order moments are 
neglected. This move closes the system but restricts the apphcability of the 
equations to behaviour on lengthscales much larger than the disk thickness. 
If the system exhibits variation on scales close to H the vertically integrated 
equations will most likely predict incorrect behaviour. 

We now have: 

DtN = -Nd^u^, (9) 
NDtUa = -Nd^i^p + - 2iVe„,^no«7 - ^f^Pap, (10) 

■^t-Pij — PijO^U^ Pj^drytli Pi^O^Uj QfSizk^O-Pkj 

— 2ejzk^oPki + Qij, (11) 

where Dt = dt + u^d^, Greek indices run only from x to y and upper case 
denotes vertical integration. Optical thickness is subsequently defined as r = 
na'^N. To complete the set, expressions for $p and $d are required. We 
approximate the central body as perfectly spherical, which accounts for the 
former potential; the latter must be obtained from Poisson's equation: 

VVi? = 4:7imGn, (12) 

in which wc may approximate n by N5{z), if it is assumed the disk is very 
thin. An additional equation for the mean vertical displacement of the disk, Z, 
can be supplied from the z component of the equation of motion (see SS85), 
and the instances of in (11) are identified with its rate of change, D^Z. But 
we shall not be investigating the vertical warping of the disk and shall assume 
symmetry about the plane z — 0. Hence Z — P^z — Pyz — 0. 

It is often illuminating to work with the internal energy and viscous stress 
equations in place of (11). The former equation may be procured by taking half 
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the trace of (11) and defining 3NC^ = Pa, where C is the vertically averaged 
velocity dispersion. A viscous stress equation proceeds by subtracting (11) 
from 25ij/3 times the energy equation, if we define the viscous stress tensor as 




We omit these details (see SS85). 
2.3 The SS85 Model 

We specialise to the modified BGK coUision term (Bhatnagar et al. 1954) 
which Shu and Stewart introduce in SS85: 



Here lJc designates the collision frequency and // is the equilibrium distribution 
function to which the coUisional process tends to push the system. This is set 
equal to a Maxwellian. Essentially the model ensures that all particles which 
undergo a collision at a certain time will be distributed immediately afterwards 
like //. Given that on the average collisions randomize velocity dispersion, this 
is not a bad approximation, but it does downplay the particles' distribution 
immediately prior to the collision: if / is very far from a Maxwellian, the BGK 
model could overestimate the isotropizing power of collisions. 

In order to account for the inelasticity of the particle collisions Shu and Stew- 
art attribute to // the same n and u as / (in order to conserve mass and 
momentum), but set its velocity dispersion cf less than c^. They then assume 
energy equipartition which permits them to write 



where the coefficient of restitution e is here understood to be a quantity aver- 
aged over collisions. The approximation of energy equipartition will worsen the 
further / departs from Maxwellian, and may explain the discrepancy between 
the results of SS85 and GT78 in the collisionless limit. 

To complete the model Shu and Stewart present an expression for the vertically 
averaged collision frequency. They find that, for a dilute, vertically isothermal 
disk with a locally Maxwellian velocity distribution. 



where fi^ is the vertical epicyclic frequency. In this case the collision frequency 
depends solely and linearly on N, and not on C. As it stands this expressions 




(2 + 



(13) 



(14) 
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fails to capture the influence of anisotropy. But it can be improved on this 
count, somewhat, via inclusion of a factor ^^NC^ jPzz-, which comes from mod- 
ifying the scale height of the disk from cj^z to Cz/^zi where Cz is the vertical 
velocity dispersion. 

The influence of self-gravity on the vertical epicyclic oscillation may be ap- 
proximated through 

If we assume (po varies slowly with x and y, then ~ d1 in Poisson's equa- 
tion, and with n oc exp(— 2;^/(2if^)) this renders ujc^ dependence on surface 
number density like: 

K 8a^ n^NsJl^ {2V2^Gm/Hnl) N. (15) 

Hence cUc can vary as steeply as iV^/^ in a self-gravitating disk, and at equilib- 
rium self-gravity enhances the collision frequency by a factor y^l -|- 2^2/7rQ~\ 
where Q is the Toomre parameter. 



2.4 The Shu & Stewart Equilibrium Disk 



The equilibrium homogeneous state of Keplerian shear can be computed an- 
alytically if the SS85 term is included, for which 



Let oUr 



Qij = uje{NC]5ij-Pij). 
u;c{N). Then, given a uniform density A^o and shear rate Uq 



— ^flox By, the uniform viscous stress tensor and the equilibrium value of e 
proceed from (11). The averaged coefiicent of restitution e depends on C 
uniquely and, once its equilibrium value is set, Cq can be deduced from the 
average of Eq.(l). The solution is 



11 + 2(u;o/^^o)2 



—3 —uq/Qo 
-cuo/Oo 2 
1 



(16) 



and 



£0 
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ll + 2(u;oM) 



(17) 



for loq = uJc{Nq) and Eq = £(Co)- 

The form of U^y demonstrates a local relationship between shear and stress. 
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Therefore one may write down the effective kinematic viscosity, 



^ (Cl\ 2{uo/no) 

[no) ll + 2(a;oW ^ ' 

as a function of ujq. Note the imphcit dependence on ujq through Cq: the 
equihbrium velocity dispersion depends on Eq which in turn is set by uuq. In 
the hydrodynamic hmit, Wo/flo — > oo, the viscosity must keep momentum 
flux constant and thus goes to zero like {cuo/^o)'^ (Chapman and Cowling 
1970). But in the collisionless limit the viscosity approaches like ujq/Q. In 
this regime the trajectories of particles between collisions are long and subject 
to epicyclic confinement. As a consequence angular momentum transport is 
inefficient. Because the effective viscosity vanishes in both hmits there exists 
a turning point in the intermediate range when the gas is hot yet not too 
rarefied. A simple calculation shows that this turning point occurs at cuq — oj, 
where 

a) = iQo ^9(l + 2/e) + y81(l + 2/e)2 + 88, (19) 
the definition of which introduces the important parameter 

e = {dlneydlnC\. (20) 

For a piecewise power- law, this parameter depends on both p and Cq/vc, but 
only weakly for the latter: we can roughly approximate e with —p. Hydrody- 
namics would suggest that one look to the region ujq > uj to observe viscous 
instability and to a;o < ct^ for viscous overst ability. Because the system is dilute 
and the particles non-spinning there is no effective bulk viscosity. 

As Goldreich and Tremaine argue, energetic instability is assured if e > 0. 
Consider a small deviation in below its equilibrium value. When e > this 
change will generate a small decrease in e and a converse increase in coUisional 
dissipation per collision (on the average). Thence coohng will shghtly domi- 
nate viscous heating and thus amplify the disturbance, initiating a runaway 
loss of energy. A small deviation above equilibrium on the other hand causes 
a decrease in collisional cooling, which will precipitate a runaway heating. If 
e — the velocity dispersion cannot adjust either way to balance the energy 
budget, and so the system will also be unstable. For e < the velocity dis- 
persion will act to return the system to equilibrium if thermally disturbed, 
and the speed of return is proportional to the magnitude of e. Isothermality 
in our model hence corresponds to e — > —oo, as then any thermal perturba- 
tion will be 'instantly' quenched. Of course 'instant quenching' is an artefact 
of the averaging process; in reality a particulate disk may return to thermal 
equilibrium no faster than the collision time, nor could an averaged e exhibit 
a steep enough dependence on c to yield e very large and negative. 

A more general perturbation including density displacements complicates the 
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simple heuristic picture sketched above, but on long lengthscales the above 
argument should hold. 



3 Linear Stability Calculation 



3.1 General Results 



In this section we estabhsh the hnear stabihty of the moment equations. We 
begin by presenting a number of results which hold for a second-order system 
derived from a kinetic equation with an unspecified collision term. 

Consider the vertically integrated equations: 

DtN = -Nd^u^, (21) 
DtUa = -dai^p + ^d) - '^ea.pQoU^ - -^d^P^p, (22) 

DtPij — Pi-yd^Uj Pj^O^Ui PijdjU^ 2rio(Qzfe-ffej ~l~ ^jzkPki) ~l~ Qiji (23) 

in a shearing sheet for a general coUision term, Qij — Qij{N, P). 

Now suppose Eq.'s (21)-(23) admit a homogeneous steady state characterised 
by the equilibrium number density A^o; Keplerian shear Uq = {0, —3D,ox/2) 
and equilibrium pressure tensor P°. Let us perturb this solution with a small 
axisymmetric disturbance: 

N = No + N{x,t), 

u = —(3^0 x/2)ey + {^(a;, t)ex + v{x, t)ej;}, 
P = P° + P(x,i), 

for \N\ < A^o, \^{u,v)\ < Qo and |P| < |P°|, then take the solution of the 
linearized Poisson equation to be 

i 2'KGm 

$n = N 

\k\ ' 

in which k is the wavenumber of the harmonic variation of N (Binney and 
Tremaine 1987). 

Next we nondimensionalise: time according to the orbital timescale, t = t*/Qo; 
space like x = (Co/f2o)x*; surface density by Nq; velocity by Co and the 
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pressure tensor by NqCq. On dropping the stars and linearising we acquire 



f) AT 




dtu = 




dtv = 


2^ ^xPxy-i 


dtPxx — 


-SP^^d^u + AP^y + Qxx, 


dtPxy — 


-2P^yd,U - P^d^V - \P,, + 2Pyy 


^tPyy = 


-P'yyd^U - 2P^^d,V - P,y + Q yy , 


dtPzz = 


-P^d^u + 4z- 



ixy 



The perturbed coUision terms take the form 



^dN)^ [dPxxJo \dP.yJo 



and self-gravity is expressed in the parameter g defined by 
which is the inverse of the Toomre Q parameter. 

We perform a Fourier decomposition into radial modes for each perturbed 
quantity so that they are oc cxp(iA;x + st), where the growth rate is s, and k — 
2/t/A is the radial wavenumber. An algebraic eigenvalue problem proceeds for 
the growth rate: Sjj Zj = s Zi, where Zi is the eigenvector {N, u, v, 

Pxxi Pxyi Pyyi Pzz) 

and Sjj is the 7x7 matrix governing the linear evolution of the system. En- 
forcing solvability returns a seventh-order dispersion relation subject to the 
equilibrium values of Pij and the derivatives of Qif 

+ (lA;^ + JA;' + K)^^ + (LA;^ + MA;^ + H)s + (PA;^ + QA;^) 
+ g\k\ {2s^ + Rs^ + (SA;2 + T)^^ + (UA;^ + \l)s^ + (WA;^ + X)s + YA;^) = 0. 

(25) 

The terms above are partitioned into those which occur independently of self- 
gravity (coefficients A-Q) and those which arise if it is included (R-Y). The 
coefficients are real and do not depend on g. 
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3.1.1 Viscous instabilities on long wavelengths 

The analysis of (25) we initially limit to wavelike instabilities on large scales, 
< /c ^ 1, i.e. wavelengths » H. And so as to capture how the various roots 
scale with wavenumber we expand s in For s — 0{k^), 

s^-^k' + 0{\k\'), (26) 

which we identify as the potential viscous instability. Its criterion for instability 
can be determined from the sign of (Q/N), which is independent of the disk's 
self-gravity. 

For s — 0{1), Eq.(25) yields the sixth order polynomial equation, 

-s^ + As^ + Cs^ + Es^ + Hs^ + Ks + N - 0. (27) 

The quadratic (1 + s^) factors (27) on account of the two identities, 

K-E + A = 0, N-H + C + 1 = (28) 

(which may be verified by computer algebra). So, to leading order, two modes 
possess growth rates s = ±i. These epi cyclic oscillations we identify as the 
potentially overstable modes. 

Without loss of generality we let s = i + p\k\ + qk^ + 0{\k\^) and return to 
(25), equating terms of order k. This obtains 

_ (V - R) + (T - X - 2)i 

^ ~ (7 + 5C - 3H - N) + (6A - 4E + 2K)i ^' 
= -ig- 

At order k'^ we collect the contributions to q from self gravity and find them 
equal to 

(15A - 6E + K - 2V + 4R) + (-11 - IOC + 3H + X - 3T)i ^ _ 1 . ^ 
(-7 + 5C - 3H + N) + (6A - 4E + 2K)i ^ ~~2^^ 

Computer algebra is required to establish both these identities. They reveal 
that the total effect of self-gravity on the overstable modes up to k"^ for A; ^ 1 
lies only in altering the wave frequency. The real part of s (emerging at order 
k'^) is unaffected. The criterion for overst ability, in fact, depends on the sign 
of 

A = (J - P - D)(7 + 5C - 3H + N) + (G - B - M)(6A - 4E + 2K). (29) 
If A > the disk is viscously overstable. 
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These results mirror those of a viscous fluid disk in the long- wavelength limit, 
as studied by Schmit and Tscharnuter 1995, Spahn et al. 2000, and Schmidt 
et al. 2001. Similarly to a fluid disk we expect self-gravity to extend the range 
of parameters for which the viscous instability and overstability occur, when 
k takes intermediate values. Without self gravity the modes of longest wave- 
length are the most susceptible to both instabilities. But inclusion of self grav- 
ity can induce the instabilities on intermediate lengthscales, Ai < A < A2 < 00, 
as well, where both Ai and A2 are in the neighbourhood of the Jeans length. 



3.1.2 Thermal Modes 

In actual fact Eq.(28) represents the decouphng of the dynamic modes from 
the thermal modes at leading order when A; <^ 1. This is because the (1,1) 
minor of Ej^ is 



E 
T 



+ 0{k) 



(30) 



where Eij is the two-by-two 'epicyclic block' and Tij is the four-by-four 'ther- 
mal block': 



E 



^ 2 ' r-j-i 9{Qxx-i Qxyi Qyyi Qzz, ^ 

9{Pxxi Pxyi Pyyi Pzz) 



4 
4020 
0-100 




We have partitioned the latter matrix into a collisional component, Cij (the 
Jacobian), and a part incorporating the effects of rotation and shear, Hij. We 
presume these components scale like |C|/|H| ~ oJc/^i i-e. the hydrodynamic 
limit will be dominated by collisional processes, and the collisionless limit by 
the anisotropizing effects of rotation and shear. Note that, to leading order in 
A;, Eq.(30) and the thermal partition also hold for a dense gas model in which 
the terms for collisional production of Pij and coUisional flux of momentum 
are arbitrary functions of A?", P^j and the rate of strain tensor. 

The 4 modes associated with the thermal block Tij include the energy mode, 
familiar from hydrodynamics (see Schmidt et al. 2001), and three others which 
emerge from the extra degrees of freedom arising from the anisotropy. In 
the coUisioness limit, |C| <C 1, we find two modes with s — ±2i + 0{\C\), 
and two with s = 0{\C\). The former pair, oscillating at twice the epicyclic 
frequency, possess, to leading order, the eigenvectors (0, 0, 0, 4, 2i, — 1, 0) and 
(0, 0, 0, — 4, 2i, 1, 0). These correspond to a perturbation of the velocity el- 
lipsoid which is comprised of the oscillation of the horizontal principal axes 
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perturbations, and their rotation with respect to the equihbrium orientation 
angle, S. One of the s = 0{\C\) modes we presume coincides with the energy 
mode sketched at the end of Section 2.3. The other is a mode affihated with 
the relaxation of anisotropy, the features of which cannot be determined until 
we specify Qij. 

In the hydrodynamic limit the thermal modes, to leading order, depend solely 
on the coUisional dynamics, their growth rates coinciding with the eigenvalues 
of Cij. These cannot be ascertained until Qij is known; but if it is presumed 
that collisions work to destroy anisotropy then we might expect the four modes 
in this limit to correspond to the hydrodynamic energy mode plus three 're- 
laxation' modes which express the decay of the deviatoric parts of the stress. 
This assumption may then be employed to determine a plausible approxi- 
mation to the collision term, which, in fact, reproduces the generic Krook 
model proposed by Shu et al. 1985 . This is clearer if we solve for and 
Uij, rather than the pressure tensor, Pij, in which case we write ^ijZj = sZi, 
where z = {N,u,v, C^,flxx,^xy,^yy)- First we adopt a simple model where 
the anisotropic parts of the stress decay independently of each other at the 
rate Ca when in the coUisional limit, but where the thermal energy mode de- 
cays at the rate ^t- Then the revised collision block, Cij will take the simple 
form: C = diag(— Ct, — Ca, — Ca, — Ca), which suggests the formula 



Qij = I F{N, C)NC^ Sij - Gij{N, U^j), 



where {dc'iF + F)q = (t, Gkk = and Gfj = — CaH^^. If Uij is small we may 
expand Gij in Iljj which reproduces a variant of SS85 (for which (a = '^c 
and F — loj^ -|- £^)/3). But as is argued in Shu ei al. 1985 the magnitude 
of the relative motions of the particles should also impact on the rate of 
thermalisation; for example, gentle or glancing collisions randomise particle 
velocities less efficiently than violent head-on collisions. The importance of 
this detail may be approximately quantified by the magnitude of G (mediated 
by e if preferred) , which leads us to connect the rate of anisotropic relaxation 
to the velocity dispersion via the first column of Gij. The coUisional block 
will hence be no longer diagonal, and then Gij will be also a function of G. 
The formula which ensues for Qij is the generalisation of Shu et a/.'s 'generic 
Krook model', with which it agrees to first order in anisotropy. 

The vertical structure of the disk may introduce anisotropy into the coUisional 
dynamics once vertical averaging is accomplished. In this case F and Gij may 
depend on II^^ as well. 
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3.2 SS85 Model Stability Analysis 



We perturb about the steady state presented in Eq.'s (16)-(17), and then 
hnearise and non-dimensionahse analogously to the previous section but with 
C — CqC* and coc — Qqco*. The stars will be henceforth dropped. 

To compute the contribution from the perturbed collision term we expand 
both £^ and cuc in Taylor series about Cq and Nq respectively. In doing so we 
have neglected any possible variation in the collision frequency due to changes 
in the disk thickness (refer to the discussion of Section 2.3), an approximation 
which simplifies the analysis. This introduces the parameters: e (cf. Eq.(20)) 
and 

uJq = {d\nUc/dhiN)o. 

If dense effects were properly included then a measure of the influence of non- 
local effects would appear in a third parameter, R = aVL/vc- But we assume 
it negligibly small. Moreover, there being no other velocity scale, all instances 
of Vc vanish in the non-dimensionalisation. Consequently, the perturbed, lin- 
earised and non-dimensionalized collision term specific to the SS85 model is 



Qij = ujo [|K + l){s^ - 1) Sij + no 1 N + luoie^ -1 + e^e) SijC^ + ujoll. 



As in the general case we have an eigenvalue problem for the growth rate 
s, subject now to the three parameters: Uq which sets the equilibrium, and 
e and uj'q which measure its linear response. The coefficients of the resulting 
characteristic polynomial in (25) are functions of these three and are listed in 
Appendix A. 

The stability problem is more tractable than Eq.(25) threatens. In fact analytic 
stability criteria can be derived for every instability the system exhibits in the 
long (and short) wavelength limits. That said the dispersion relation could 
be simphfied substantially if we put into use Eq.(14) and the average of (1), 
in which case lOq — 1 and e ^ —0.234, —0.20, —0.19. But we leave these 
as free parameters. Firstly, doing so provides a means by which the collision 
frequency's enhancement by the disk's self-gravity can be included. This effect 
should ensure Uc is superlinear in N and thus 3/2 > cUq > 1. Secondly, leaving 
cUq and e open frees the system to exhibit the viscous overstability. 



3.2.1 The Viscous Instability 



To begin we restrict the analysis to large scales, < /c -C 1, and suppose 
s = C»(A;2). Given e < 0, it follows from (A.14) that N < 0; so, from (26), the 
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criterion for instability is Q > which, on employing (A. 16), obtains 

, -e(l + a;g)(ll + 2o;g) 

° 18a;2 + e (1 + u;i){ll - 2Lui) ' ^ ^ 

when ujq > uj, where uj is the nondimensionahsed turning point of u, introduced 
in (19). This assures us that the instability is indeed the kinetic analogue 
of the viscous instability: if we substitute the expression for our particulate 
disk's effective kinematic viscosity (Eq.(18)) into the hydrodynamic criterion 
for viscous instability, 

we recover precisely (31). Moreover we find that the growth rate is 

s = -(Q/N) + 0(\k\^) = -Suo {(3 + l)k'' + 0(1^;^ (32) 

on using (A. 16) and (A. 14) again, which is exactly the growth rate of the 
viscous mode as computed by hydrodynamics. 



Two curves of marginal viscous stability in the (a;o,ci;Q) plane are plotted in 
Fig. 2. They diverge as ujq approaches uj from above, and asymptote to a;Q = 1 
as a;o — > oo- A simple calculation demonstrates that for e < —9/11 a marginal 
curve will asymptote to 1 from above, and for e > —9/11 from below. In the 
hydrodynamic regime (for a dilute ring) we may make the identification uj'q — 
—(3 + 0{uJq'^). Therefore the kinetic criterion of viscous instability recovers the 
familiar hydrodynamic criterion, (5 < —1. 

When e — > 0, both the curves' turning point and uj drift to the origin, rendering 
more of parameter space unstable. At e = the entire quadrant is unstable, 
which coincides with the energy instability sketched at the end of Section 2.4. 
On the other hand as e — > — oo (the isothermal limit) we find uj — > a 
finite value. 

Our discussion of uJcS form at the end of Section 2.3 would suggest that 
appropriate values of uj'q lie between 1 and 3/2. Thus e > —9/11 will ensure 
that such disks are viscously unstable for all collision frequencies above a 
critical value close to uj. A quick calculation shows this value to be (— lle/(9 + 
lle))^/^. The data offered by experiment suggest e < —9/11, and thus indicate 
a dilute disk composed of hard, icy boulders falls into this category. 

As the assumption of diluteness fails in Saturn's rings (as discussed earlier), 
the validity of these stability predictions is unclear. In fact, theories and simu- 
lations of 'dense' and low-optical thickness equilibria suggest that the viscous 
instability should not appear (Araki and Tremaine 1986 and Wisdom and 
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Fig. 2. Marginal Curves for the viscous instability in the (ooq, ijj'q) plane for <C 1 and 
different e representing the two regimes of e > —9/11 and e < —9/11. The regions 
above the curves are unstable. The dotted lines encompass the region 1 < loq < 3/2 
in which we expect to find equilibria characteristic of a dilute ring. 

Tremaine 1988). Only in the simulations of Lukkari 1981 and Salo 2001 has 
the instability been observed, though the former used a scaling of the elasticity 
law which did not account for nonlocal effects, and the latter was for a 2D 
disk. So no physically realistic 3D simulation to date has exhibited the insta- 
bility. This we presume is because either the parameters are not appropriate 
{aD,/vc or T too big) or the lengthscales of the computations have not been 
sufficiently large to capture the unstable modes. 

A final observation: if we let uJq = 1 and solve for e we produce the stability 
criterion: e > — (9/ll)a;Q/(l +u!o), which is qualitatively the same as that of 
Stewart et al. 1984 who solve a simplified kinetic model for only the optical 
thickness and thermal energy. So it appears that reduction is justified here. 



3.2.2 The Viscous Overstability 

In the long wavelength limit, A = describes the curves of marginal viscous 
overstability. Substituting the required coefficients obtains 



2(c/i +92e + g3e^) 



(33) 



crit 



9(l + a;2)(/ii + /i2e + /i3 e^) 
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where the g^'s and /I's are the polynomials: 



gi{uJo) = 9(-50886 - 75014o;^ - 2AAA0uj^ - U79lu^ + 139u^ + 10cj(5°), 
c,2(o;o) = 6(990 + 5855a;^ + 9991a;^ + 60Uu^ + 959ojg + 73uj^^ + 2u;^^), 
gsM = ^1(1 + a;o)'(-1851 - ISllo;^ - lOOa;^ + 4a;^), 
hiiuo) = 9(-2178 + 833ujI + 891a;^ + 156a;^ + 8uj^), 
h2{uJo) = a;o(705 + 396a;g - 351a;^ - 42u;^), 
/laK) = 2a;o(l + ^o)'(-63 + llc^o + 2c^o)- 

The modes are unstable when uj'q > u![.J.^^ , and the denominator in (33) negative 
(corresponding to u!o> uj); and when cu'q < ui'^rit, and the denominator positive 
{cuq < a;). Thus the stability curve comprises two branches in the {cuqjUJq) 
plane, as plotted in Fig. 3. 

The right branch approaches the line o^q = (| + |e~^) from below in the 
hydrodynamic limit, a;o — > oo, and only falls into the first quadrant of the 
{ujo,u!q) plane if e < —3. This is a value experiments show to be too low for 
hard ice particles. Moreover the maximum cJq for which overstability in the 
hydrodynamic limit is possible is 2/9. Thus a disk also needs to exhibit a 
strongly sublinear dependence of lJc on N for overstability to develop, which 
is implausible. However, a number of formal results follow by making the 
identification /3 = —ouq; we then obtain the hydrodynamic criterion (3>(3c — 
— (| + |e~^); and if, in addition, the isothermal limit, e — oo, is enforced, 
we recover precisely the overstability criterion of Schmit and Tscharnuter for 
an isothermal fluid disk without bulk viscosity: (5 > —2/9. 

The left branch is associated with the monotonically increasing dependence 
of equilibrium viscosity with when cjq < cD- Fig. 3 shows these curves to 
diverge when uj ^ uj, and intersect the a;o = axis at a finite value of a;o, which 
some algebra reveals to be 514/99 — (20/297) e. Subsequently the minimum 
uj'q for which left-branch viscous overstabihty is possible is ^ 5.192 — a highly 
superlinear dependence of Uc on N. It follows that for all ujq and e there exists 
a 'channel', 2/9 < cJq < 5, in which overstability is forbidden. The equilibria 
of a dilute ring lie well within this region, hence viscous overstability will not 
occur in a dilute ring. 

3.2.3 Thermal Modes 

We complete the A; <^ 1 analysis by examining the four modes embodied in the 
SS85 thermal block T, analogous to that appearing in Eq.(30). We find that 

the cooling component of this matrix is the expected, C = diag(— r(72 , —ujq, —ujq 
where we have defined the collisional cooling rate to be F = ^C'^NuJc{1 — e^), 
and the subscript indicates partial differentiation with respect to squared 
velocity dispersion. 
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Fig. 3. The two branches of the marginal curves of viscous overstability in the 
(a;o,Wo) plane for different e and with A; <C 1. The left branch intersects the loq = 

514/99 — (20/297)e, and the right branch asymptotes to the line 

< 3/2. 



axis at a;' 
, ,/ 



- 

2,2 
9 ^ 3' 



'^). Probable values of lOq lie in the dotted channel 1 < ojj 



In the hydro dynamic limit, cjq — oo, the heating contribution to T is sub- 
dominant and the stress components, Hxx, ^xy and Hyy possess (to leading 
order) a decay rate equal to uq (representing the 'relaxation modes'), while 
the velocity dispersion, C^, decays at the rate Tc-^ (corresponding to the mode 
in which the energetic stability of the equilibrium is expressed) . 

In the converse collisionless regime, a;o ^ 1; we recover the complex conjugate 
pair we met earlier. Their growth rates are 



s — 



±2i - ^(100 + e)u;o + 0{u;l)] + 0{\k\), 



and so are stable for all e > —100, an upper bound which is much larger than 
the values of e we deem feasible. These two modes are analogous to those of the 
viscous overstablity, only becoming unstable when the thermal restoring forces 
are sufficiently strong (i.e. e negative and large). The mechanism animating 
them we presume works to quench perturbations of the horizontal axes of the 
velocity ellipsoid, but may be so effective that the resulting 'overshoot' may 
instigate an oscillation of growing amplitude. The instability can therefore be 
thought of as an 'anisotropic overstability'. 

For general ujq there is a single 'relaxation mode' with growth rate s — —cuq + 
0{\k\), composed exclusively of U^x and Uyy perturbations in the ratio 4 to 
1. This mode deforms the velocity ellipsoid, cither vertically squashing or 
vertically dilating it, but while maintaining C constant. Thus it is a pure 
rearrangement of anisotropy. Factoring out this mode from the characteristic 
polynomial of Tij leaves a cubic in s with a single real solution corresponding 



22 



to the growth rate of the energy mode. This is negative if: 



|u;o(ll + 2u;^) 



2uJn 



11 + 2uj^ 



> 0. 



(34) 



The square bracketed term we identify as the partial derivative of the energy 
balance with respect to squared velocity dispersion, C^, because it is composed 
of the derivative of the coUisional cooling function, F, and the viscous heating, 
|z/o. Thus thermal stability on large scales depends solely on the response of 
the energy equilibrium to small fluctuations in velocity dispersion. When we 
substitute Tq^ into Eq.(34) the condition simplifies to e < 0, which justifies 
the heuristic argument offered at the end of Section 2.4. 



3.2.4 Instabilities on Intermediate Lengthscale 

Now we solve for general k. Caution should be exercised in interpreting the 
results of this analysis, considering the omission of the third order moments. 
These will certainly dominate short-scales but should also play some part 
on intermediate lengths. We will discuss their effects in more depth in the 
hydrodynamic analysis of Section 4. 

Let us begin without self-gravity, i.e. let = 0. We seek the marginal stability 
of non-oscillatory modes, and hence let s = 0. The curves of marginal stability 
for a mode of wavenumber k are determined by 

PA;^ + Q = 0, 

which derives from Eq (25). The criterion Q = describes the marginal stabil- 
ity of the longest waves and coincides with the criterion of viscous instability. 
The criterion P = corresponds to marginal stability on the shortest wave- 
lengths, and is associated with an energy mode. Because we assume e < 
thermal stability will hold on long wavelengths (see the analysis in Section 
3.2.3) but on shorter scales this need not be true. Outside the long wave- 
length regime density perturbations can no longer be factored out and the 
energy balance's dependence on N plays a part. We shall discover in the hy- 
drodynamic analysis that this energy mode attempts to keep both pressure 
and viscous stress constant, and its stability is related subsequently to Field's 
criterion (Field 1965) . That said the omitted terms representing the heat flux 
we presume enter at higher order than k^ and should extinguish this instability 
on sufficiently short scales. On intermediate scales instabihty is guaranteed if 
Q/P < 0, which refers to a hybrid mode of viscous and thermal effects. Again, 
omission of the third order moments will probably affect its stability. 

With self-gravity added the curves of marginal stability for a non-oscillatory 
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mode are described by 

Pk^ + Yg\k\ + Q = 0. 

Solvability for \k\ furnishes the criterion > g'^^.^^ = 4QP/Y^, or equivalently, 
Q < Qcrit, where Q is the Toomre parameter. The critical value g^crit is a 
complicated function of luq, lUq and e, which we list in (B.l). 

In Fig. 4 curves of marginal stability are plotted in {u!q,ujq) space for = 1 
and e = —0.2. Three curves are plotted. Equilibria in the region above the 
dotted graph will develop viscous instability on all wavelengths above a crit- 
ical lengthscale. This is precisely the same curve we encountered in the long 
wavelength analysis, as the most unstable wavelengths are the longest. The 
solid graph represents marginal viscous/gravitational instability, above which 
equilibria become unstable. The unstable region includes that of the viscous 
instability, which is a special case. Modes which are viscous/gravitationally un- 
stable possess wavelengths in an interval (Ai, A2) bracketing the Jeans length; 
modes which are viscously unstable possess A2 = 00. The marginal curve of 
viscous/gravitational instability asymptotes to a line ujq — l—g^/3 for ujq large, 
a fact easily verified from Eq.(B.l). It follows that in the hydrodynamic limit 
the instability criterion can be written as f3 < (3^ = ~(1 ""fi'^/S), which agrees 
with Schmit and Tscharnuter's criterion (Eq.(35) in Schmit and Tscharnuter 
1995) for an isothermal disk. This fact suggests that the thermal properties of 
the ring play a neghgible role on the stability criterion of the most unstable 
mode in the collision-dominated, hydrodynamic limit. 

The dotted-dashed curve in Fig. 4 represents the quasi-thermal instability, 
which affects all modes with wavelengths below a critical intermediate value. 
Perhaps we can then characterise it, like the viscous instability, as a special 
case of gravitational/viscous instability but in which Ai = 0. Inclusion of a 
heat flux term would undoubtedly prevent Ai being zero. So in reality scales 
sufficiently small will be stable. 

Suppose that we take ou'q — 3/2, e = —0.234 and g = 1. For k small the SS85 
model predicts that viscous instability will be present in those parts of a dilute 
ring in which there are more than 3.5 collisions per orbit (i.e. ujq > 0.555), a 
critical limit which corresponds to an optical thickness of 0.218 if we employ 
expression (14). If, however, we include cUcS self-gravity factor, introduced in 
equation (15), the viscous instability can occur for lower optical thicknesses; 
for example if g = 1 instability occurs for r > 0.089. Hence only the most 
tenuous of dilute rings will be stable. 

The SS85 model also suggests that on intermediate scales quasi-thermal or 
viscous/gravitational instability will occur in parts of the ring in which there 
are less than about 33 collisions per orbit (tuo ^ 5.2318), which corresponds to 
a critical r of approximately 2 (or less if a;c's self-gravity factor is included). 
This last result is dubious, of course, because of the omission of the third 
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order moments. The prevalence of quasi-thermal and viscous/gravitational 
instability on intermediate scales is a question left open until the next section. 

We will not undertake an analogous analysis of the viscous overstability in 
the intermediate range, it being algebraically intensive and not especially il- 
luminating. In summary, the inclusion of self gravity narrows the 'channel 
of stability' in the (cjq, cJq) plane, but not sufficiently to endanger the stabil- 
ity of equilibria with appropriate values of uJq. Instability on those regions 
of parameter space rendered 'newly' unstable by self-gravity disturb modes 
on intermediate scales only; thus the 'old' boundaries of stability circumscribe 
regions in which viscous overstability operates on arbitrarily long wavelengths. 

It should be mentioned that in the k ^ 1 regime the system not only evinces 
the unphysical extension of the quasi-thermal instability but also two pairs of 
unstable sound waves, which we also regard as unphysical. That said these only 
appear for parameter values incongruent with a dilute ring of icy particles: the 
first pair of modes may only be unstable when e < —65, and the second when 
uj'q < 0.98 for e and ujq sufficiently small. Thus the thermal instability will be 
the primary problem if one sought to undertake a non-linear computation of 
the time-dependent equations. 



4 Comparison with Viscous Fluid Disk 



In this section we compare the linear analysis of our kinetic system with that 
of hydrodynamics. The motivation for doing so is to distinguish the effect of 
the viscous stress equation, which should not only capture anisotropic effects 
but also the non-Newtonian behaviour of the stress in the range of low and 
intermediate ujq. We produce an analysis analogous to Schmidt et al. 2001, and 
thus solve the mass conservation equation, the Navier-Stokes equation, the 
temperature equation and Poisson's equation for a self-gravitating, viscous 
fluid disk in a shearing sheet. We omit however Schmidt et al.'s nonlocal 
contributions to pressure and viscosity. So as to best mimic our kinetic disk 
the transport coefficients have been set to those computed by the SS85 kinetic 
equilibrium. Hence the coefficient of bulk viscosity and heat diffusivity, n, are 
zero for the moment, the latter corresponding to the omission of the third 
order moments. The shear viscosity is 

llQ2 + 2u;,(7V)2' ^^^> 

where T = is temperature. In addition we incorporate the 'cooling func- 
tion', r, into the temperature equation which shall account for coUisional 
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Fig. 4. Marginal Curves for the viscous instability, quasi-thermal instability and 
viscous/gravitational instability for general k in the (wcWq) plane for e = —0.2 
and g = 1, as predicted by the SS85 kinetic model. Regions above the curves are 
unstable. Thus equilibria in Region A are stable. Region B circumscribes equilibria 
susceptible to gravitational/ viscous instability on modes with k G (fci, /c2)- Equihb- 
ria in Region C feature viscous instability affecting modes with A; < /ca. In Region 
D both quasi-thermal instability and viscous instability are apparent, and so every 
mode is unstable. And in Region E both quasi-thermal and gravitational/viscous 
instability exist on all wavenumbers greater than a critical value k4. The curve 
of viscous instability as3niiptotes to = 1 for large wq, while the curve of vis- 
cous/gravitational instability goes to cjq = 1 — 

energy losses, 

T = ^uj,NT{l-e{Tf). (36) 
The fluid equations in the shearing sheet read: 

DtN = -NdkUk, 
N DtUi = -N dii^p + ^d) - diP + dkTlik, 
iNDtT^-PdkUk + ^-r, 

where pressure is defined by P = NT and the rate of viscous heating is 
^ = UkiSki, where Sij — eij — dkUkSij/3 designates the rate of deformation 
tensor. The stress is Newtonian, so 11^/ = 2Ni/Ski, and is determined from 
Poisson's equation (12). It is understood these equations, and the fields which 
appear in them, are vertically averaged. 

The stationary state is represented by: the equilibrium surface number density 

A^'o, the Keplerian shear u = —^VIqX By, and the equilibrium temperature 
To, which is computed from the energy balance. We perturb, linearise and 
non-dimensionalise this state as earlier, expanding both v and F in Taylor 
series. On assuming wavelike perturbations oc e**"^'*^^ the system yields the 
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non-dimensional fourth-order dispersion relation: 

+ {Et + li^k^} 

+ {1 - 2g\k\ + (I + fE^uo + 3z/T lyo)^ + \vl k^}s' 

+ {E^{\ - 2g\k\) + (Et - z/o + 2z/T + O^o - frN] 

- 2^z/o \kf + d^T^^o + 7z/tz^o + %)vQk^} s 
+ {2(1 + OFt - 2rNZ/T/z^o - 2g (fFx + |z/t) 

+ (KFt - Tn) - 1(1 + Oi^o + I^^t) k^j^o k' = 0, (37) 

where ^ = (^In z//ain iV)o, = {du/dT)o, and = i(rT - *t) which 
is 2/3 times the partial derivative of the system's net rate of cooling with 

respect to temperature. Note that ,^ is a partial derivative and thus differs 
from (3 introduced in Section 1. If heat diffusivity k is included other terms 
materialise, the ones of relevance to our analysis being: 

i/Q k'^ilnk^ - l9K\k\^ + 2(C + l)^^;^). 
For more details see Schmidt et al. 2001. 

In order to complete the linear theory of the fluid analogue we replace the 
four extra parameters appearing in the dispersion relation with 

_ 9^0(1 + 0^^) 1 9-2e(l + a>g) 

" 2(ll + 2o;g)' - ll + 2a;g ' ^^^^ 



°(ll + 2a;2)' ^ ll + 2u;2 ^ J 



4-1 Modes on Large Lengthscales 



In the long-wavelength limit, A; <^ 1, and assuming s — 0{k'^), the viscous 
mode has growth rate: 

,^_^il±^Ei^^^l2^k2^o{k') 

E'-£ 

On applying Eq.'s (38)-(39), curves of marginal viscous instability are derived 
which are identical to those predicted by the kinetic model in Eq.(31). Also if 
we let s — 0{1) we find the long- wavelength thermal mode: 

s^-ET + 0{k'^). 

which reproduces Eq.(34). 
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The growth rate of an overstable mode Schmidt et al. show to be 
s^i- ig\k\ + |(7i - + |(7^ + |z/o + 3Cz/o)/c^ 

where 

1 



7i = 



7r = 



1 + 



1 



2 1 + 3^ 2 

3 2 " 3 



— 3i^T I N -C'T ^ i^O -C'T — g 



in the absence of nonlocal terms. Hence the curves of marginal overstability 
proceed from 

_ -4rN - 4:Et + 7u + AE^Uq + 3z/t(4 - 4£^TrN - ^EtUq) 
~^ ^ 9uo{^ + 2E^ + 3EtUt) ■ 

which translates, in the parameters of our dilute ring, to 



crit 



(ll + 2a;g)(/i + e/2 + e^ h 



(40) 



where 



h = 45(11 + 2cj, 
Z5 = -72a;o'(l + ^o), 



2\2 
2\2 



k = 6(121 + 183a;^ + 66lu^ + 4lu^), 
/4 = 9(-1331 + 132cj^ + 16lu^), 

Z6 = 4a;2(-ll + 2c^o)(l + ^o)'- 



These curves resemble roughly those of the kinetic model, consisting of two 
branches: one tightly bunched near the collisionless limit and one which ex- 
tends over the remainder of parameter space (See Fig. 5). In the limit ujq ^ oo 
the graphs asymptote, like the kinetic curves, to lu'q = {^ + |e~^). However for 
intermediate and lower values of ujq the two accounts disagree markedly. The 
kinetic and fluid graphs diverge at different values of the collision frequency, 
ujq. And the left branch intersects the a;o = axis at lJq ^ —0.556 — 0.074e. 
Hence for e > —1.874 overstability occurs for all uj'q to the left of the sin- 
gularity. This contradicts the kinetic model which predicts that no viscous 
overstability can occur for 2/9 < lu'q < 5.192. Thus a dilute ring of icy parti- 
cles with e = —0.234, lo'q — 3/2, and g — 1 manifests viscous overstability for 
all ujQ < 2.2022, or r < 0.541 

This is an important qualitative difference, and which we argue results from 
the non-Newtonian nature of the viscous stress, in particular its nonlocality 
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Fig. 5. Marginal curve of viscous overstability predicted by hydrodynamics for 
e = —0.2 and k <^ 1. The right branch asymptotes to the hne u'q = 2/9 + 2/(3e), 
for large ujq (in agreement with the kinetic model) but diverges to +oo for an inter- 
mediate value of loq, as opposed to the kinetic model which diverges to -co. The 
left branch differs markedly, intersecting the loq = axis for negative Uq and thus 
ensuring a large sector of parameter space is unstable. The channel 1 < cJq < 3/2 
indicates the region in which we expect appropriate equilibria to fall. 

in time. This effect interferes with the phase synchronisation of the stress and 
density oscillations, and thus the mechanism sustaining the overstable modes. 
It is less clear how anisotropy is at work here, but we conclude it plays a 
smaller effect in the linear theory. This effect we expect to be less appreciable 
for a dense gas, as the coUisional contribution to the viscous stress, which is 
local in time, will be of order or greater than the free-streaming component, 
which is non-local in time. 



4-2 General Stability Criterion for Non- Oscillatory Modes 



We derive a relation describing the marginal stability of non-oscillatory modes 
from (37) if we set s = 0. It follows that instability occurs on an interval of k 
if the following equation admits two solutions for \kc\, 



i(rT-rN)-|(i + 0^o 



|Z/T 



kl ~ 2g 



3^ T 



+ [2(i + e)rT-2rNZ/TM] = o, (41) 



Eq.(41) is the counterpart of (32), the coefficients of k"^, \kc\ and A;^ corre- 
sponding to P, Y and Q, respectively. Thus the coefficient of determines 
the onset of viscous instability and dominates in the long wavelength limit. 
The coefficient of \kc\ introduces self-gravitation. The coefficient of k'^ mani- 
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fests thermal effects, and is in fact equal to the derivative of 2/3 times the net 
rate of cooling, = |(r — with respect to temperature, if both pressure 
and the viscous stress are kept constant. To show this we must include the 
variations induced by temperature and density displacements upon the shear 
rate, A, which compensates for the variation in v induced by the same per- 
turbations. If Hxy is to be kept constant, given a fluctuation in N and T at 
constant P, then 

f)A 

^ = iKM-(i + 0], 

and therefore 




= §(rT-rN)-i(i + 0^^o + i^T. 



It follows that the criterion for quasi-thermal instability on very short scales 
(when Eq.(41) is dominated by the thermal term) is a modification of Field's 
criterion (dE/dT)p < (Field 1965). But on intermediate scales, the poten- 
tially unstable mode will be an amalgam of effects associated with the thermal 
balance, self-gravitation and angular momentum transport, all constrained to 
maintain constant pressure and stress. 

The marginal curves which follow from Eq. (41) match those of SS85 (described 
in Fig. 4); the viscous/gravitational graph nearly coincides, and the quasi- 
thermal graph is qualitatively the same. Here the discrepancies which do exist 
between the two models we attribute to the breakdown of the assumptions of 
both on intermediate scales. Their details we omit. 

Heat diffusivity introduces the additional terms — |(7k|A';|"' + ^K,k^ which con- 
tribute at sufficiently small scales. Including heat conduction will also incur 
a k"^ term, 2(^ -|- l)K.k^, representing the coupling of thermal transport with 
rotation and viscosity. Hence the criterion of intermediate scale instabihty will 
be altered significantly. This may be quantified by a critical heat diffusivity, 
Kcrit below which quasi-thermal instability can emerge on a band of k, the 
upper limit depending sensitively on k,. An expression for this quantity in the 
non-self-gravitating case is presented in Eq.(C.l). We find this critical value to 
be very small indeed; for instance, when e = —0.234, lJq — 1 and (Vo < 0.7 then 
'^crit ^ 0.07. The dilute kinetic theory asserts that u k, (Cowhng and Chap- 
man 1970), so for (35) only dilute rings of large optical thicknesses (r > 3) 
should feature the quasi-thermal instability. 

When the destabihsing effect of self-gravity is added, however, KciH can ap- 
proach feasible values. Instability in this case we name viscous/gravitational, 
though in fact the mode in question is fueled by viscous, gravitational and 
thermal processes. In order to establish how this instability impacts on the 
various scales we plot curves of marginal stability in the {u!o,u'q) plane for 
given e, k, and k, as described by Eq.(C.2) (Fig. 6). When k < 1 we find that 
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Fig. 6. Marginal curves of non-oscillatory instability for three representative values 
of k as predicted by the hydrodynamical model, when e = —0.2, g = 1, and k = I. 
Equilibria situated above the A; ~ curve are unstable on the longest scales, but 
those above the A: = 3 curve are stable for scales fc > 3. Marginal curves for larger k 
disappear into the top right. The area outside the region bounded between the two 
branches of the dotted-dashed curve are stable on scales k > 0.5 

the marginal curves fall into three categories corresponding to long, interme- 
diate and short scales. For scales k < kic we obtain curves resembling those of 
viscous instability (see Fig. 2). But for scales kic < k < k2c there appear two 
branches, the right one familiar from Fig. 4, the left arising from non-zero k. 
Instability emerges on equilibria situated between the two. We conclude that 
heat diffusion confines instability on intermediate scales to a smaller area of 
parameter space than predicted otherwise; in particular, equilibria with low 
ujo are stabilised. For short waves, k > k2c, the marginal curves mimic in 
shape those of long wavelength, but as k gets large 0;^^,;^ > k^. Therefore if 
small scales are to be rendered unstable the system needs to exhibit an in- 
appropriately large value of cu'q. As expected, heat diffusivity stabilizes short 
scales. Expressions for kic ^-nd ^2c be obtained by examining the real 
roots shared by the numerator and denominator of (C.2). For n = 1, they 



are approximately equal to g± ^g'^ — (6/1 l)a;, which lie on either side of half 
the Jeans length. The fluid model then predicts that a dilute ring, with self 
gravity {uJq = 3/2) and k = 1, is stable on scales A* < 27ry^2/3, or about 5 
times the disk semi-thickness. 

It is likely a kinetic theory, accounting for heat conductivity adequately, would 
predict similar results, though this cannot be checked until such a model is 
devised. However, hydrodynamics shows clearly that if we wish to study be- 
haviour on intermediate and short scales, it is necessary to include in some 
way the third-order moments. 
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5 Compcirison with the Goldreich &c Tremaine Kinetic Model 



This section outlines the hnear stabihty of the GT78 model, so as to com- 
pare the modified Krook collision term with that of a formal Boltzmann-type 
account. 

Goldreich and Tremaine take the distribution function, /, to be a three- 
dimensional Gaussian in velocity space with its axes corresponding to the 
principal axes of the pressure tensor, pij. Assuming symmetry about the plane 
z — 0, the principal axis system is denoted by the unit vectors ei, e2 and 
^3 — ^z, so that Bz sin 5 — Bx x ei, where 5 denotes the orientation angle. The 
principal values ofpij are accordingly defined by pkk = ncf., where k = 1, 2, 3, 
and the summation convention is not applied. The velocity distribution is then 



This form results in a much simplified collision term; see Eq (33) in GT78. 
However to facilitate better the comparison with SS85, we assume the value of 
s, which appears in the coUisional kinematics, is a quantity averaged over col- 
lisions, and not a constant. Also we note that the implicit collision frequency 
in the resulting Qij is oc fir Ck/c^ and thus is dependent on the velocity disper- 
sions, being 'weighted' differently for each direction. It does not capture self 
gravity or filling factor effects in its present form- these requiring a superlinear 
dependence on r. 



5.1 Equilibrium Solution 



For the equilibrium analysis we solve the four components of the pressure 
tensor equation (10) in the principal axis frame, cf. Eqs (37) in GT78. Given 
T, these determine e, 5, and the axis ratios of the velocity ellipsoid, C2/C1 and 
C3/C1. 

The results of our numerical root finding reproduce those of GT78 and can be 
observed in Fig. 7, in which we plot the e — r law only. The SS85 relationship is 
plotted for comparison. Eq.(14) has been used to relate r and Uc- As mentioned 
earlier the discrepancy at low r is marked. For more discussion of this see SS85 
and Shu et al. 1985. 




(42) 
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Fig. 7. The energy equilibrium relations between the coefficient of restitution, eq, 
and optical thickness, ro, as determined by the SS85 and GT78 models. 

5.2 The GT78 Linear Stability 

In order to model the action of self gravity on the collision frequency we 
re-parametrise the disk equilibrium and replace optical thickness by a;o, as 
expressed in (14). For the purposes of the linear stabihty analysis we then 
treat Uc as a supcrlincar function of N. This we admit is somewhat crude 
but should approximate to an acceptable degree self- gravity's enhancement of 
the collision frequency, and certainly sharpen the comparison with the SS85 
model. 

We now perturb about the steady state calculated in the previous section. Wc 
then linearise and non-dimensionalise, though on account of the complicated 
form of the collision term this process is rather involved: 

(1) First, we must calculate, in the principal axis frame, the three collision 
term perturbations, Qkk, in terms of perturbations of r, S, and the prin- 
cipal velocity dispersions, Ck- 

(2) Then the Ck and 6 perturbations are converted to those of the Pij com- 
ponents using coordinate transformation rules. 

(3) Lastly Qij is written in terms of Qkk using the transformation rules. 

In this scheme s is assumed a function of C^, where = l{cl + 0^ + c|), and 
is Taylor expanded about the steady state. After it is assumed that the per- 
turbations are wavelike and axisymmetric, a seventh-order dispersion relation 
emerges of the same structure as (25) and dependent on the three parameters, 
cuq, cu'q and e, alongside k. 

Given e, the marginal curves of the viscous instability and overstability are 
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Fig. 8. Comparison between the SS85 and GT78 kinetic models with respect to the 
marginal curves of viscous instability, for e = —0.2 and e = —1.5. Regions above 
the curves are unstable 

numerically determined in the {cuq, cu'q) plane. They are plotted in Figs 8 and 
9. 

The curves for the viscous instability share the same quahtative features as 
those of the SS85 model. Though the former deviate somewhat in shape, their 
singularities approximately coincide, they asymptote to the line cJq = 1, and 
the critical value of e (which determines whether a curve approaches this line 
from above or below) is roughly —9/11. This suggests that not only is the 
viscous mode resistant to the effects of non-locality in time (as proven in the 
previous section) but also to the details of the collision term, in particular to 
anisotropy, for which the Goldreich and Tremaine model accounts somewhat 
better. 

Both branches of the curves of overstability agree in their general features 
with those of SS85. In the large ujq limit, the right branches asymptote to 
Uq = est; in the isothermal limit this constant is 2/9, according with the SS85 
and fluid models. However the critical e, below which marginal curves can 
cross into the first quadrant, is closer to —2.4 than —3. Also the left branches 
intersect the uq = axis at slightly lower values: stability is assured in the 
region left of the singularity if uj'q < 4.674, as opposed to ^ 5.192. We conclude 
that the linear stability of the viscous overstability is a little more sensitive 
to the details of the collision term, and hence anisotropy — though not nearly 
as much as to non-locality in time. As previously mentioned, the mechanism 
of overstability relies on the synchronisation of two oscillations, that of the 
acoustic-inertial wave, and that of a forcing which issues from the accompa- 
nying variation in angular momentum flux. The instability is vulnerable to 
disruption in the phase and in the magnitude of this forcing. A stress free 
to oscillate independently certainly impacts on the former. Anisotropy may 
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Fig. 9. Comparison between the SS85 and GT78 models with respect to the marginal 
curves of viscous over stability when e = —0.2. Regions above the left branch are 
unstable while regions below the right branch are unstable. 

impact on the latter. 

Lastly we resolve the potential 'anisotropic overstable' modes for low luq and 
find them unstable for e < —48, far larger than predicted by the SS85 model 
(though still unfeasible) . This emphasises again the fact that the two models 
differ most in their treatment of anisotropy. We conclude that the less accurate 
SS85 model, only linear in Hij, fails to fully model the deformations of the 
velocity ellipsoid. 



6 Conclusion 



In summary, we have calculated the linear stability of a patch of a dilute 
planetary ring in a shearing sheet employing a second-order kinetic theory. 
We have defined 'diluteness' as the insignificance of non-local effects, whose 
importance we can roughly quantify by aQ/c for lower optical depths. But 
throughout Saturn's rings the velocity dispersion is low enough and particle 
size large enough, for this quantity to be non-negligible, thus we recognise 
that the applicability of a dilute model is restricted. Nevertheless this work 
reveals the interesting effects anisotropy and non-Newtonian stress have on the 
viscous instability and overstability. It also provides a preliminary framework 
in which a linear analysis of a dense gas may be undertaken. 

We adopted a system of vertically integrated moment equations derived from 
Shu and Stewart's inelastic generalisation of the BGK equation in their 1985 
paper. The ensuing system manifests both the viscous instability and over- 
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stability in the long- wavelength limit, and their extensions induced by self- 
gravity on intervals of intermediate wavelength. It predicts viscous or vis- 
cous/gravitational instability in dilute rings where particles collide about 3.5 
times or more per orbit (optical thicknesses greater than 0.2, or less depending 
on the importance of self-gravity on the collision frequency) . Viscous oversta- 
bility does not feature in dilute regions. 

These results were compared with a hydrodynamic model incorporating the 
kinetic theory's effective viscosity (evaluated at equilibrium) and cooling func- 
tion. We find that the criterion for viscous instability is precisely the same as 
that for the kinetic model, but the two criteria diverge in the intermediate 
range of collision frequency for the viscous overstability. According to the ki- 
netics, overstability is suppressed in this regime, a result in contradiction with 
the fluid model which predicts the instability to be widespread. We conclude 
that hydrodynamics presents a misleading impression of the viscous oversta- 
bility's prevalence on small and intermediate ranges of cJc/f^Oi foi' ^ dilute ring. 
This arises from the increasingly non-Newtonian nature of the stress in this 
regime, an effect the kinetic account captures but which the hydrodynamical 
cannot. We stress that this incongruence will be most marked for a dilute ring; 
as the coUisional stress tensor is local in time, we assume the denser the ring 
the better the hydrodynamic approximation. 

We also undertook a linear stability analysis of the Goldreich and Tremaine 
1978 kinetic model, employing a Boltzmann-type collision term with a triaxial 
Gaussian, which should factor in anisotropic effects more thoroughly. The 
criteria for viscous instability and overstability mirror qualitatively that of the 
SS85 model, and thus suggest that the linear dynamics of a disk are insensitive 
to the the precise details of the collision term. In particular, the criterion of 
viscous instability matches closely that of the SS85 model. The mechanism 
of this instability would thus seem relatively robust, not only resistant to the 
effects of non-locality in time but also of anisotropy. The discrepancy between 
the marginal curves of viscous overstability computed by the SS85 and GT78 
models is a little more pronounced. We conclude the viscous overstability is 
relatively fragile for small and intermediate tUc/^o — very susceptible to non- 
locality in time and slightly altered by anisotropy 

On intermediate and short scales the SS85 model also predicted the emergence 
of a quasi-thermal instability which the hydrodynamical analysis revealed to 
be associated with a mode orienting the thermal and density perturbations in 
such a way to keep the equation of motion undisturbed. This mechanism thus 
mimicks that operating in Field's thermal instabihty (Field 1965). We find that 
the instability should occur on optical thicknesses less than 2. However hydro- 
dynamics shows that this instability is severely sensitive to heat conduction; 
K need not be large to extinguish it entirely. Thus the omission of third-order 
moments in the kinetic model bestows on it a prominence that it does not de- 
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serve. When self-gravity is included intermediate scale instability can occur for 
much larger k, which in this case we refer to as a 'viscous/gravitational insta- 
bility'. Self-gravity, in addition to its important dynamic role on intermediate 
scales, will enhance the collision frequency, and render the dependence of cUc 
on N as steeply as N^^"^. In general this will make the system more susceptible 
to both the viscous instability and overstabihty on long lengthscales. 



Two possibly fruitful avenues which proceed from this work are: linear stability 
calculations of dense and/or spinning particle disks, and nonlinear studies of 
the evolution of the viscous instability and overstabihty using the second-order 
moment equations. 



In principle the techniques employed here can be set upon kinetic models 
of dense and spinning rings. Detailed formulations have been developed in 
Shukhman 1984, Araki and Tremaine 1986, Araki 1988 and Araki 1991, but, 
unmodified, their algebraical truculence surely prohibits all but equilibrium 
calculations. The success of the SS85 model encourages us that a simphfication 
of the dense disk collision terms could provide a model that captures the 
essential physics, and also expedites linear and nonlinear analysis. The results 
obtained with such a model would be directly pertinent to N-body simulations 
and Saturn's B-ring. It could also ascertain the validity of the fluid description 
in the dense regime. 



It would be instructive also to determine the impact of the nonlinear evo- 
lution of the second-order moments on the development of the instabilities 
discussed, as compared to the Newtonian fluid model, which so far has only 
been explored in the isothermal case (Schmit and Tscharnuter 1999). They 
may elicit interesting nonlinear behaviour from the viscous overstabihty. But 
an understanding of the small-scale instabilities, even if unphysical, is essential 
for the sucessful numerical implementation of these equations. 
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A Appendix: Dispersion Relation Coefficients 



The coefficients appearing in the dispersion relation of the SS85 model (see 
Eq.(25)) are: 



A 
B 
C 
D 

E 
F 
G 



2 e(l+a;g)-9(6 + a;g) 
36 



11 + 2tug' 
2eujl{l + ul) - (55 + 52c;2 ^ g^4-) 



9(ll + 2a;2)' 



3(ll + 2cjg' 
12(10 + cc;2)' 



11 + 2cjg 

20e(10 + IIujI + u;^) + 3c^^,(198 + 2,QujI) 

-3(4084 + 1037a;^ + 58u;^) 
e(13 + IQo;^ + 6a;^) - 3(80 + 31a;^ + 2a;^) 



2\2 



(11 + 2ujI) 
2 

9(ll + 2a;g) 

, .2 



e Uq (553 + 599u;o + mul) - 27 e c^q (1 + c^g) 



+4a;^ (2349 + 4056,;^) + (-3960 - 13785a;^ - 2856a;^ - 132t,;^ 



H = 



1+^0 



3(ll + 2a;2) 



ec;^(17 + 2ujI) - (132 + 18t,;^^ 



(A.l) 
(A.2) 
(A.3) 

(A.4) 
(A.5) 

(A.6) 



(A.7) 
(A.8) 
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I 



2e (700 + 849a;^ + 156a;^ + lu^) 



9(11 + 2iji)^ 

+3u;'q (1980 + 558a;^ ^ 35^4^) _ 3(21460 + 73446^^ ^ 739^4 ^ ggcu^) 



(A.9) 



J 



2ci;o 



e (1 + a;o)(80 + 32u;o + ^^0(478 - 99a;^)) 



9(11 + 2cu2)2 

-3(10u;^ + a;^(1977 - 819a;^) + u;^(303 - 1350;^) + 4(250 + 99u;o)) 



K 



3(11 + 2u;^) 

9(11 + 2ljo)^ 



e (11 + 17a;^ + 60;^) - 3(44 + 25uj^ + 2uj^) 
e (1 + cu^)(2240 + 20cu^ + cj^(451 - 36uj'q) - 270o;^, 



(A. 10) 
(A.ll) 



-6(5156 + 5uj^ + u;^(1578 - 459a;o) - 2034cUo - 9a;^(-17 + 2c;;)) 



M 



3(11 + 2uiy 



e (1 + a;o)(25 + 2uj^ + u;^(45 - 30^;^) + 33^;^) 



+6(cu^(-4 + 60;^,) - 2(8 + 334) + + 48^)) 



(A.12) 

(A.13) 
(A.14) 



4u;g 



3(11 + 2u;i)^ 



e (1 + a;^)(176 + 2uj^ + a;^(43 - 12a;^,) - IUuJq) 



24 

(11 + 2a;o)2) 

4a;o 
3(ll + 2cx;§) 

S = 2 + ^^, 
11 + 2cug' 

4(1 + u^i) 



e{l + u',){2u;l( 



+72u;'q{25 + iu;^) 
ll(l + u;(,)}-18u;^,^ 



(A.15) 



e(l+a;^) -9(6 + 0;^) 



T = 
U = 
V = 



11 + 2uji 

3(11 + 2cuiy 
2a;o 
3(ll + 2u;g 



e a;Q — 22 — SuJq 



-2e (10 + llu;^ + cu^) + 3(316 + 77u;^ + 4uj^) 
e (11 + 17a;2 + q^4^ _ 3(^44 + 250;^ + 2uj^) 



(A.16) 
(A.17) 

(A.18) 

(A.19) 

(A.20) 

(A.21) 
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2\2 



3(11 + 2cj^ 
3(11 + 2ujIY 



e (49 + 53a;^ ^ 4^4^, _ g^^gg ^ 260;^ 



e(29 + 31u;^ + 2ul) - 18(10 + 



(A.22) 
(A.23) 
(A.24) 



B Appendix: Critical Self Gravity Constant 



We present the expression for the critical value of g above which gravita- 
tional/viscous instability arises in the SS85 kinetic model. See Section 3.1.1. 



a;, 



5'crit 



'I [hi{uJo, e) + h2{uo, e) cv'q + ^3('^o, e) 



121 



(11 + 2a;g)3(e(29 + Slcu^ + 2uj^) - 18(10 + cu^)) 



(B.l) 



where 



hi = -e'(16 + a;^)(ll + 13a;^ + 2a;^)^ 

h2 = e(ll + 130;^ + 2uj^) [e(-62 - 29a;^ + 35uj^ 



h. = -6 



216a;^(25 + 4a;^) - 6e(-550 - 481cj2 + gicj, 



- 18(100 + 32u;2 ^ 
22cj^) 



+6^(1 + 6^2)2 (_209 + IGuul + 4cu^) 



C Appendix: Critical Coefficient of Heat Diffusivity 



Here we hst the critical k below which quasi-thermal instability occurs on a 
band of wavelengths for a viscous fluid model without self-gravity. 

_ 9(1 + QVg + (1 + Ot/o(4rN + 4rT - 9t/T) - sr^uT 

12(1 + 0^^0 

J(4rN + 9(e + iHM + 1)^0 - ^T)((e + i)rT^o - Tn^t) 

12(1 + 0% 

(C.l) 



When self-gravity is present, the marginal curves in the {ujq, uj'q) plane for given 
k, K and e may be determined from Eq.(41) with heat flux terms added. For 
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a dilute ring with an equilibrium computed by SS85, we obtain 



11 + 2ul 



3 [IScug + 33a;ofc^ + c^oe(ll + 9a;^ - 2uj^) + nk'^^Au^ - 121)] 



X 



a;o(l + a;^)(3 + k^)e + 2gkuJo{9 - e(l + a;^)) - Kk'^{ll + 2a;^)(3 - 2gk + fc^ 

(C.2) 
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